/************************************************************************
*************************************************************************
*************************************************************************

gnr_jpe.do

This program runs the GNR procedure using their own code. I've made one
major change: the user can now specify whether to assume a linear or nonlinear
specification for the Markov process of productivity. Note that the program
requires the data to be xtset before running. It also creates several variables
while running, and will throw an error if those variables already exist.
Finally, it produces the elasticities of output as variables.

9 April 2015: 	Created
28 July 2015: 	Modified to use syntax
20 July 2016:	Added a version that uses only lagged labor as an instrument.
27 Feb  2017:	When the linear option is specified, the program now stores the
				estimated AR coefficient on productivity


Arguments:
varlist			- 	Level of output,labor,capital,flexible inputs, share of inputs (in that order)

Optional:


The program creates new variables that contain the estimated elasticities
on each of the inputs.


*************************************************************************
*************************************************************************
************************************************************************/



/************************************************************************
Main program
************************************************************************/

cap program drop gnr_jpe
program gnr_jpe
	syntax varlist, [ linear laginst]
	
	
	
	local word = 0
	foreach loc in yg_level l_level k_level i_level si_level {
		local word = `word' + 1
		local `loc': word `word' of `varlist'
	}
	
	gen yg = ln(`yg_level')
	gen l = ln(`l_level')
	gen k = ln(`k_level')
	gen i = ln(`i_level')
	gen si = ln(`si_level')

	gen ll=l*l
	gen kk=k*k
	gen ii=i*i
	gen lk=l*k
	gen li=l*i
	gen ki=k*i
	gen lki=l*k*i
	*gen lll=ll*l
	*gen kkk=kk*k
	*gen iii=ii*i


	/************************************************************/
	/**  The code below runs the share regression using a      **/ 
	/**  a (log)polynomial approximation.  This step     	   **/
	/**  recovers the output elasticity of flexible inputs 	   **/
	/**                                                        **/
	/**  The initial values are set using an OLS regression    **/
	/**  of shares (si) on a polynomial in inputs.  If the     **/
	/**  non-linear least squares procedure (nl) fails to      **/
	/**  converge, the initial values can be changed.	   **/
	/************************************************************/

	// Obtain initial values from OLS refression of flexible inputs shares on all productive inputs //
	regress si l k i ll kk ii lk li ki lki if si~=. & l~=. & k~=. & i~=.
	matrix test = e(b)
	predict temp_c if si~=. & l~=. & k~=. & i~=.
	replace temp_c = temp_c - _b[_cons]

	//Correct constant term for initial values //
	egen mtemp_c = min(temp_c)
	scalar ntemp_c=mtemp_c
	drop temp_c mtemp_c
	scalar ntemp_c=-ntemp_c + 0.1

	// Using initial values from OLS, regress shares using log polynomial approximation //
	#delimit;
	nl ( si = ln({g0=ntemp_c} + {gl=_b[l]}*l + {gk=_b[k]}*k + {gi=_b[i]}*i + {gll=_b[ll]}*ll + 
	   {glk=_b[lk]}*lk + {gli=_b[li]}*li + {gkk=_b[kk]}*kk + {gki=_b[ki]}*ki +
	   {gii=_b[ii]}*ii + {glki=_b[lki]}*lki) ) if si~=. & l~=. & k~=. & i~=., iter(100);
	#delimit cr

	// Construct predicted elasticities of output with respect to intermediate inputs //
	predict elasm_gnr if l~=. & k~=. & si~=. & i~=.

	// Perform mean adjustment of predicted elasticities //
	predict eg if l~=. & k~=. & si~=. & i~=., resid
	replace eg=-eg
	egen mexp_eg=mean(exp(eg))
	replace elasm_gnr=elasm_gnr-ln(mexp_eg)
	replace elasm_gnr=exp(elasm_gnr)

	// Rename and save mean adjusted coefficients of estimation //
	mat beta=e(b)
	svmat double beta
	ren beta1 g0
	ren beta2 gl
	ren beta3 gk
	ren beta4 gi
	ren beta5 gll
	ren beta6 glk
	ren beta7 gli
	ren beta8 gkk
	ren beta9 gki
	ren beta10 gii
	ren beta11 glki
	foreach var of varlist g0-glki {
		egen s`var'=mean(`var')
		drop `var'
		ren s`var' `var'
		replace `var' = `var' / mexp_eg
	}	
	clear matrix

	// Compute the integral //
	gen integ_G_I = g0+gl*l+gk*k+gll*ll+gkk*kk+glk*lk + (gi*i+gli*li+gki*ki+glki*lki)/2 + gii*ii/3
	replace integ_G_I=integ_G_I*i 

	// Construct script_Y //
	gen vg = yg - eg - integ_G_I

	// Generate lagged variables //
	//tset id time
	gen vg_1=L.vg
	gen l_1=L.l
	gen k_1=L.k
	gen ll_1=L.ll
	gen kk_1=L.kk
	gen lk_1=L.lk
	gen vg_2=L2.vg
	gen l_2=L2.l
	gen k_2=L2.k
	gen ll_2=L2.ll
	gen kk_2=L2.kk
	gen lk_2=L2.lk

	/************************************************************/
	/**  Now to recover the remaining coefficients associated  **/
	/**  with capital and labor (the constant of the PDE)      **/
	/**                                                        **/
	/**  The initial values are set using an OLS regression    **/
	/**  of vg on log capital and loglabor                     **/
	/************************************************************/

	// Set initial values using OLS estimation of script_Y on capital, labor, their squares and interaction term //
	reg vg l k ll kk lk 
	matrix test2 = e(b)
	matrix test3 = test2[1,1..5]
	matrix test3[1,1] = test2[1,"l"]
	matrix test3[1,2] = test2[1,"k"]
	matrix test3[1,3] = test2[1,"ll"]
	matrix test3[1,4] = test2[1,"kk"]
	matrix test3[1,5] = test2[1,"lk"]
	matrix drop test2

	// Recover capital and labor coefficients using gmm procedure //
	// gmm_prod function computes moment conditions using productivity innovation and capital and labor inputs //
	
	//If we want to impose the AR assumption, do the linear program; otherwise, do the usual program
	if !mi("`linear'") {
		gmm gmm_prod_lin`laginst' if vg!=. & l!=. & k!=. & vg_1!=. & l_1!=. & k_1!=., one nequations(5) parameters(al ak all akk alk) from(test3) winitial(identity) rhs(vg l k ll kk lk vg_1 l_1 k_1 ll_1 kk_1 lk_1) conv_maxiter(100)
	}
	else {
		gmm gmm_prod`laginst' if vg!=. & l!=. & k!=. & vg_1!=. & l_1!=. & k_1!=., one nequations(5) parameters(al ak all akk alk) from(test3) winitial(identity) rhs(vg l k ll kk lk vg_1 l_1 k_1 ll_1 kk_1 lk_1) conv_maxiter(100)
	}


	mat beta=e(b)
	svmat double beta
	ren beta1 al
	ren beta2 ak
	ren beta3 all
	ren beta4 akk
	ren beta5 alk
	foreach var of varlist al-alk {
		egen s`var'=mean(`var')
		drop `var'
		ren s`var' `var'
	}	

	/************************************************************/
	/**  Generate productivity in logs and levels, and compute **/
	/**  the output elasticities of labor and capital.         **/
	/************************************************************/

	// generate log productivity //
	//gen logomega=vg-al*l-ak*k-all*ll-akk*kk-alk*lk
	//gen omega=exp(logomega)
	
	//If the AR option is specified, generate the AR estimate
	if !mi("`linear'") {
		
		tempvar omega omega_1
		gen `omega' = vg-al*l-ak*k-all*ll-akk*kk-alk*lk
		gen `omega_1' = l.`omega'
		reg `omega' `omega_1'
		gen ar_gnr 		= _b[`omega_1']
	}
	else {
		gen ar_gnr 		= .
	}

	// Generate output elasticities of labor and capital inputs //
	gen elasl_gnr=gl*i+2*gll*li+glk*ki + gli*ii/2 + glki*ki*i/2 + al + 2*all*l + alk*k
	gen elask_gnr=gk*i+2*gkk*ki+glk*li + gki*ii/2 + glki*li*i/2 + ak + 2*akk*k + alk*l
	
	order elas*, alpha last
	drop yg-alk
end




/************************************************************************
The function evaluator program for the second stage
This version uses a third-order polynomial for the Markov process
************************************************************************/
cap program drop gmm_prod
program gmm_prod
	version 11

	syntax varlist [if], at(name) rhs(varlist)
	
	local m1: word 1 of `varlist'
	local m2: word 2 of `varlist'
	local m3: word 3 of `varlist'
	local m4: word 4 of `varlist'
	local m5: word 5 of `varlist'

	local vg: word 1 of `rhs'
	local l: word 2 of `rhs'
	local k: word 3 of `rhs'
	local l2: word 4 of `rhs'
	local k2: word 5 of `rhs'
	local lk: word 6 of `rhs'
	local vg_1: word 7 of `rhs'
	local l_1: word 8 of `rhs'
	local k_1: word 9 of `rhs'
	local l2_1: word 10 of `rhs'
	local k2_1: word 11 of `rhs'
	local lk_1: word 12 of `rhs'

	tempname al ak al2 ak2 alk
	scalar `al'=`at'[1,1] 
	scalar `ak'=`at'[1,2]
	scalar `al2'=`at'[1,3]
	scalar `ak2'=`at'[1,4]
	scalar `alk'=`at'[1,5]

	tempvar w w_1 w2_1 w3_1 csi

	quietly gen double `w'=`vg'-`al'*`l'-`ak'*`k'-`al2'*`l2'-`ak2'*`k2'-`alk'*`lk' `if'
	quietly gen double `w_1'=`vg_1'-`al'*`l_1'-`ak'*`k_1'-`al2'*`l2_1'-`ak2'*`k2_1'-`alk'*`lk_1' `if'
	quietly gen double `w2_1'=`w_1'*`w_1' `if'
	quietly gen double `w3_1'=`w2_1'*`w_1' `if'
	quietly reg `w' `w_1' `w2_1' `w3_1' `if'
	quietly predict `csi' `if', resid
	
	quietly replace `m1'=`l'*`csi' `if'	
	quietly replace `m2'=`k'*`csi' `if'	
	quietly replace `m3'=`l2'*`csi' `if'	
	quietly replace `m4'=`k2'*`csi' `if'	
	quietly replace `m5'=`lk'*`csi' `if'	

end


/************************************************************************
The function evaluator program for the second stage
This version uses a linear (AR) polynomial for the Markov process
************************************************************************/
cap program drop gmm_prod_lin
program gmm_prod_lin
	version 11

	syntax varlist [if], at(name) rhs(varlist)
	
	local m1: word 1 of `varlist'
	local m2: word 2 of `varlist'
	local m3: word 3 of `varlist'
	local m4: word 4 of `varlist'
	local m5: word 5 of `varlist'

	local vg: word 1 of `rhs'
	local l: word 2 of `rhs'
	local k: word 3 of `rhs'
	local l2: word 4 of `rhs'
	local k2: word 5 of `rhs'
	local lk: word 6 of `rhs'
	local vg_1: word 7 of `rhs'
	local l_1: word 8 of `rhs'
	local k_1: word 9 of `rhs'
	local l2_1: word 10 of `rhs'
	local k2_1: word 11 of `rhs'
	local lk_1: word 12 of `rhs'

	tempname al ak al2 ak2 alk
	scalar `al'=`at'[1,1] 
	scalar `ak'=`at'[1,2]
	scalar `al2'=`at'[1,3]
	scalar `ak2'=`at'[1,4]
	scalar `alk'=`at'[1,5]

	tempvar w w_1 w2_1 w3_1 csi

	quietly gen double `w'=`vg'-`al'*`l'-`ak'*`k'-`al2'*`l2'-`ak2'*`k2'-`alk'*`lk' `if'
	quietly gen double `w_1'=`vg_1'-`al'*`l_1'-`ak'*`k_1'-`al2'*`l2_1'-`ak2'*`k2_1'-`alk'*`lk_1' `if'
	quietly reg `w' `w_1' `if'
	quietly predict `csi' `if', resid
	
	quietly replace `m1'=`l'*`csi' `if'	
	quietly replace `m2'=`k'*`csi' `if'	
	quietly replace `m3'=`l2'*`csi' `if'	
	quietly replace `m4'=`k2'*`csi' `if'	
	quietly replace `m5'=`lk'*`csi' `if'	

end






/************************************************************************
The function evaluator program for the second stage
This version uses a third-order polynomial for the Markov process

This one uses only lags of labor as instruments (no longer
imposing that labor is completely determined the year
before).
************************************************************************/
cap program drop gmm_prodlaginst
program gmm_prodlaginst
	version 11

	syntax varlist [if], at(name) rhs(varlist)
	
	local m1: word 1 of `varlist'
	local m2: word 2 of `varlist'
	local m3: word 3 of `varlist'
	local m4: word 4 of `varlist'
	local m5: word 5 of `varlist'

	local vg: word 1 of `rhs'
	local l: word 2 of `rhs'
	local k: word 3 of `rhs'
	local l2: word 4 of `rhs'
	local k2: word 5 of `rhs'
	local lk: word 6 of `rhs'
	local vg_1: word 7 of `rhs'
	local l_1: word 8 of `rhs'
	local k_1: word 9 of `rhs'
	local l2_1: word 10 of `rhs'
	local k2_1: word 11 of `rhs'
	local lk_1: word 12 of `rhs'

	tempname al ak al2 ak2 alk
	scalar `al'=`at'[1,1] 
	scalar `ak'=`at'[1,2]
	scalar `al2'=`at'[1,3]
	scalar `ak2'=`at'[1,4]
	scalar `alk'=`at'[1,5]

	tempvar w w_1 w2_1 w3_1 csi

	quietly gen double `w'=`vg'-`al'*`l'-`ak'*`k'-`al2'*`l2'-`ak2'*`k2'-`alk'*`lk' `if'
	quietly gen double `w_1'=`vg_1'-`al'*`l_1'-`ak'*`k_1'-`al2'*`l2_1'-`ak2'*`k2_1'-`alk'*`lk_1' `if'
	quietly gen double `w2_1'=`w_1'*`w_1' `if'
	quietly gen double `w3_1'=`w2_1'*`w_1' `if'
	quietly reg `w' `w_1' `w2_1' `w3_1' `if'
	quietly predict `csi' `if', resid
	
	
	quietly replace `m1'=`l_1'*`csi' `if'	
	quietly replace `m2'=`k_1'*`csi' `if'	
	quietly replace `m3'=`l2_1'*`csi' `if'	
	quietly replace `m4'=`k2_1'*`csi' `if'	
	quietly replace `m5'=`l_1'*`k_1'*`csi' `if'

end


/************************************************************************
The function evaluator program for the second state
This version uses a linear (AR) polynomial for the Markov process

This one uses only lags of labor as instruments (no longer
imposing that labor is completely determined the year
before).
************************************************************************/
cap program drop gmm_prod_linlaginst
program gmm_prod_linlaginst
	version 11

	syntax varlist [if], at(name) rhs(varlist)
	
	local m1: word 1 of `varlist'
	local m2: word 2 of `varlist'
	local m3: word 3 of `varlist'
	local m4: word 4 of `varlist'
	local m5: word 5 of `varlist'

	local vg: word 1 of `rhs'
	local l: word 2 of `rhs'
	local k: word 3 of `rhs'
	local l2: word 4 of `rhs'
	local k2: word 5 of `rhs'
	local lk: word 6 of `rhs'
	local vg_1: word 7 of `rhs'
	local l_1: word 8 of `rhs'
	local k_1: word 9 of `rhs'
	local l2_1: word 10 of `rhs'
	local k2_1: word 11 of `rhs'
	local lk_1: word 12 of `rhs'

	tempname al ak al2 ak2 alk
	scalar `al'=`at'[1,1] 
	scalar `ak'=`at'[1,2]
	scalar `al2'=`at'[1,3]
	scalar `ak2'=`at'[1,4]
	scalar `alk'=`at'[1,5]

	tempvar w w_1 w2_1 w3_1 csi

	quietly gen double `w'=`vg'-`al'*`l'-`ak'*`k'-`al2'*`l2'-`ak2'*`k2'-`alk'*`lk' `if'
	quietly gen double `w_1'=`vg_1'-`al'*`l_1'-`ak'*`k_1'-`al2'*`l2_1'-`ak2'*`k2_1'-`alk'*`lk_1' `if'
	quietly reg `w' `w_1' `if'
	quietly predict `csi' `if', resid
	
	quietly replace `m1'=`l_1'*`csi' `if'	
	quietly replace `m2'=`k_1'*`csi' `if'	
	quietly replace `m3'=`l2_1'*`csi' `if'	
	quietly replace `m4'=`k2_1'*`csi' `if'	
	quietly replace `m5'=`l_1'*`k_1'*`csi' `if'

end
